* Gary LaMotte's Mathematica solutions to exercises from the 2012 Gravitational Wave Astronomy Summer School

*DAY 5  PART3.*

*a.*

samprate=16384;   nyquist = samprate/2;

dataC =  Flatten[Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day5dataC.txt", "Table"]];

dataD =  Flatten[Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day5dataD.txt", "Table"]];

tdataC =  Table[(r-1.0)/samprate  , {r,Length[dataC]}];tdataD =  Table[(r-1.0)/samprate  , {r,Length[dataD]}];
ndataC = Length[dataC];dataCsecs = ndataC/samprate;  ndataD = Length[dataD];
dataDsecs = ndataD/samprate;dataCsecs = ndataC/samprate;

*b.*

crosscorr = dataC.dataD

crosscorr = dataC.dataD\[Conjugate]

shiftcorr =Table[dataC.RotateLeft[dataD,i*8192],{i,1,10}];

rmsshiftcorr=RootMeanSquare[shiftcorr]

SNR = crosscorr/rmsshiftcorr

*c.*

*FOR DATA C:  Apply Square PassBand Filter Subprogram:*

Clear[lowcorn,upcorn, filterdataC]

filterdataC[lowcorn_,upcorn_]:=
(start = {1,1}; pass = PadRight[start, ((upcorn - lowcorn)*dataCsecs)+1, 1]; rtside = PadRight[pass, ndataC-((lowcorn*dataCsecs)-1),0]; sqfilt = PadLeft[rtside,ndataC];  
sinc = InverseFourier[sqfilt]; realsinc =  2*Re[sinc];ListConvolve[realsinc/Sqrt[Length[realsinc]],dataC,1])

lowcorn =  40;

upcorn = 200;

filtdataC = filterdataC[lowcorn,upcorn];

*ALSO,FOR DATA D:  Apply Square PassBand Filter Subprogram:*

Clear[lowcorn,upcorn, filterdataD]

filterdataD[lowcorn_,upcorn_]:=
(start = {1,1}; pass = PadRight[start, ((upcorn - lowcorn)*dataDsecs)+1, 1]; rtside = PadRight[pass, ndataD-((lowcorn*dataDsecs)-1),0]; sqfilt = PadLeft[rtside,ndataD];  
sinc = InverseFourier[sqfilt]; realsinc =  2*Re[sinc]; 
ListConvolve[realsinc/Sqrt[Length[realsinc]],dataD,1])

lowcorn =  40;

upcorn = 200;

filtdataD = filterdataD[lowcorn,upcorn];

filtered`cross`corr = filtdataC.filtdataD

shiftcorr =Table[filtdataC.RotateLeft[filtdataD,i*8192],{i,1,10}];

filtered`rms`shiftcorr=RootMeanSquare[shiftcorr]

SNR = filtered`cross`corr/filtered`rms`shiftcorr

*d.*

*Non-Optimal Part*

fftC = Fourier[dataC];   fftD = Fourier[dataD];

f`cross`corr = Re[fftC.fftD\[Conjugate]]

shiftedD = Table[RotateLeft[dataD,i*8192],{i,1,10}];

shiftcorr = Re[fftC.Transpose[(Fourier[shiftedD])\[Conjugate]]];

f`rms`shiftcorr=RootMeanSquare[shiftcorr]

nonoptimal`freq`SNR = f`cross`corr /f`rms`shiftcorr

*Optimal Part*

noise =  Flatten[Import["C:\\Users\\Gary LaMotte\\Documents\\Shawhan GW Lectures and Exercises\\day3asd.txt", "Table"]];

noisespectrum = noise[[1;;262144]];

noisepsd = Abs[(noisespectrum )^2];

halffreqarr =  Table[(arr/dataDsecs)^-3  , {arr,nyquist*dataDsecs}];

reversehalffreqarr = Reverse[halffreqarr];

joinedfreqarr = Join[halffreqarr ,reversehalffreqarr];

*The following two inputs restrict freqarr values to have a maximum of 1*

max[x_]:= If[x>1, 1, x]

freqarr = max/@joinedfreqarr;

kernal = freqarr*(1/noisepsd);

optf`cross`corr = Re[fftC.(fftD\[Conjugate]*kernal)]

optf`shiftcorr = Re[fftC.(Transpose[(Fourier[shiftedD])\[Conjugate]]*kernal)];

optf`rms`shiftcorr=RootMeanSquare[optf`shiftcorr]

optimal`freq`SNR = optf`cross`corr /optf`rms`shiftcorr